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We report on a strong coupling approach (on-site Coulomb repulsion, U larger than the nearest- 
neighbour hopping energy \t\) to the Hubbard model. Starting from the Hubbard operators which 
diagonalize the interaction term, we generate a hierarchy of composite operators from the equations 
of motion. Using the Hubbard operators as a basis, we are able to compute the associated Green 
functions including the anomalous Green functions which describe pair formation. We show explicitly 
that these anomalous Green functions are non-zero in the A x 2_ y 2 channel; however, the entities that 
pair up are not single electron-like particles but rather composite excitations (which we call cexons) 
made out of an electron and a hole on nearest-neighbour sites. Cexons are fermionic in nature as they 
have spin 1/2 and also have unit charge. Our calculations of the chemical potential reveal that negative 
compressibility in the 2D Hubbard model and composite excitation pairing are intimately connected, 
namely, the larger the negative compressibility, the larger the pairing amplitude. Our observation 
of negative compressibility in the under-doped regime is consistent with phase segregation or stripe 
formation in the normal state. While pairing ameliorates the negative compresssibility, it does not 
eliminate it entirely. In addition, we find that the anomalous correlation functions are particle-hole 
symmetric and exhibit a maximum at a doping level of roughly 10% as measured from half-filling. For 
U — 8\t\, the onset temperature for pair formation is . 02 1 £ | . The effect of nearest-neighbour Coulomb 
repulsions is discussed. 



I. INTRODUCTION 

After a protracted quest 0] for superconductivity in 
the Hubbard model, it is surprising that the question 
is still open. This state of affairs has arisen primarily 
because the parameter space in which superconductivity 
is anticipated to reside is precisely the strong coupling 
limit, U 3> \t\, in which traditional perturbative schemes 
fail. As a result, progress on this question has relied pre- 
dominantly on exact diagonalization or Quantum Monte 
Carlo (QMC) simulations on finite systems. While early 
work [jjj showed promising results on the possible onset 
of superconducting pair correlations in the 2D Hubbard 
model, the most recent numerical work || shows that 
true superconductivity with at least algebraic decay of 
anomalous correlations is absent. In light of such re- 
sults, it has also been suggested that some combination of 
phonons || and or next nearest-neighbour hopping pj are 
needed for superconductivity to survive. Further, Kuroki 
and Aoki Q have suggested that QMC simulations have 
yet to access the energy scale of order O.Olt where strong- 
coupling superconductivity || is expected to occur. In 
addition, Beenen and Edwards J?J have analyzed the 2D 
Hubbard model with an equations of motion technique 
considerably enhanced by Mancini, Matsumoto, and co- 
workers ||] and shown that pairs with d x i_ y i symmetry 
emerge. However, this work has key weaknesses: the on- 
set temperature and the window of doping over which 
d x 2_ y 2 pairing occurs are highly sensitive to the decou- 
pling scheme of the anomalous Green functions. For ex- 



ample, at least order of magnitude variations in the onset 
temperature were observed for U = 8\t\ and a de- 
coupling scheme which is supposedly exact in the limit 
of U — > oo resulted in a near vanishing of pairing for 
U > 10\t\. 

Despite these difficulties, the Hubbard model remains 
central to such strongly correlated problems as super- 
conductivity in the cuprates as well as the organic con- 
ductors. In fact, in so far as the motion of holes in the 
CuC>2 planes in the cuprates can be modeled realistically 
H by the Hubbard model in the vicinity of half-filling, 
the microscopic underpinnings of the pairing mechanism 
in these materials should arise from this model. Because 
D=2 is the marginal dimension for superconductivity, it 
might be that while the pairing mechanism originates 
within a single layer, interlayer tunneling is needed to 
maintain long-range phase coherence. Hence, the ques- 
tions that face the Hubbard model are three-fold: 1) 
does local pairing of the right sort obtain to describe 
the cuprates, 2) what entities are involved in the pairing, 
and 3) is long-range phase coherence present? In this 
work we develop a strong-coupling approach that is suf- 
ficient to answer the first two questions. Our approach 
is based on a simple observation first made by Hubbard 
|iof . Namely, that the bare electron annihilation opera- 
tor can be written as a sum of two composite operators 

C-ia Cier^-i — er T Cier(l Hi — o) ^]ia H~ da- (1) 

The quantities described by r\i a and £j CT are composite 
excitations. The r\ excitation describes an electron re- 
stricted to move on sites already occupied with an elec- 
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tron of opposite spin whereas £ demands that there be no 
prior occupancy on the site. A key feature these opera- 
tors possess is that they exactly solve the t — Hubbard 
model. Consequently, an equation of motion approach 
based on these operators will lead naturally to an ex- 
pansion in t/U. Such an expansion is ideally suited for 
the strong coupling regime U \t\ in which both the 
cuprates as well as the organic conductors reside. It is 
this approach that we pursue here. 

When the kinetic energy is treated as a perturbation, 
hybridization is introduced into the Hubbard atomic or- 
bital basis. As a result, £j and r\i are no longer the eigen- 
operators. The new type of excitations can be thought 
of heuristically as resonant valence-bond hybrids. To de- 
termine what operators create such excitations, it is suf- 
ficient to construct the equations of motion for the Hub- 
bard operators. The relevant operators have the kinetic 
energy, t as a coefficient. As we will see, the simplest 
terms that arise are of the form, r\ij a = a a nj T with (ij) 
nearest neighbours. This operator is the generator of a 
composite excitation (or cexon) that is restricted to live 
on sites i with a nearest neighbour site j occupied. Such 
composites have unit charge as can be seen directly from 
the commutator, 

[Vija-,nitr>] = 8<r,<T>5iir)ij<r. (2) 

Also, quite obviously, its spin is 1/2. In fact, all the 
composite excitations generated in our scheme have unit 
charge. Hence, spin and charge are coupled. However, 
while the rjij excitation is fermionic in character, it is 
not confined to a single site. Rather, it is an exci- 
tation that lives on neighbouring sites. In traditional 
perturbative schemes, such as Fermi liquid theory, simi- 
lar composite operators naturally appear. However, in 
such approaches, the operators that are generated do 
not describe new excitations but rather dress the non- 
interacting quasiparticles. In this approach, the new op- 
erators do in fact describe fundamentally new excitations 
flTf as they give rise to new peaks in the spectral func- 
tion. The appearence of new peaks in the spectral func- 
tion is incompatible with a Fermi liquid description of 
the composite excitations. 

It might be argued that the emergence of cexons in 
this model and the unrecoverability of Fermi liquid the- 
ory are inevitable in an approach based fundamentally 
on the atomic limit of the Hubbard model. However, any 
theory of high T c must be based on a starting point that 
captures the essence of strong-correlation physics. While 
the atomic limit might seem extreme, it is justified in this 
context as U ^S> t in high T c systems. Consequently, it 
stands to reason that the physically-relevant excitations 
should bear a greater family resemblance to the atomic 
limit eigentates than they would to the plane-wave exci- 
tations which populate the non- interacting limit. Cexons 
are the simplest local excitations that arise once the Hub- 
bard operators are hybridized on nearest-neighbour sites. 



To determine if such excitations really do constitute new 
particles in the strongly-coupled limit, dynamical correc- 
tions must be included to the theory we present below. 

As in earlier composite operator approaches |7],||,[l2|,|l3| , 
we find that the traditional correlation function, (c^Cji), 
does not determine the pairing gap in the 2D Hubbard 
model. Rather 

% = (c iT c a n ir ) = (rhjlVijl) (3) 

is the relevant anomalous correlation function that de- 
termines pairing in the d-wave channel. The emergence 
of this correlation function as the d-wave order parame- 
ter in the Hubbard model and the seemingly innocuous 
rewriting (overlooked previously in similar Green func- 
tion analyses in terms of the r\ij a operators is par- 
ticularly illuminating. Pairing based on the r\ij a excita- 
tions is the new feature which we develop explicitly in 
this work. This rewriting signifies that composite exci- 
tations rather than electron-like particles, as would be 
the case in a traditional Fermi liquid, are involved in the 
pairing process in a strongly-correlated system with re- 
pulsive interactions. In addition, the fact that the new 
composite operators describe the pairing lends further 
credence to their correspondence with new excitations in 
the strongly-correlated regime. While it is well known 
that quasiparticles in strongly-correlated systems do not 
resemble well electron-like excitations, as shown, for ex- 
ample, by the absence of a sharp peak at the Fermi sur- 
face in angle-resolved photoemission experiments Jl4[ on 
the under-doped normal state of the cuprates, the tra- 
ditional mechanism invoked to explain this fact is spin- 
charge separation |],|l5| . In 1-dimension |ll| spin-charge 
separation is on firm footing. In addition, even in a sin- 
glet BCS superconductor, spin-charge separation occurs 
naturally as the Cooper pairs carry charge but no spin, 
whereas the quasiparticles are spin 1 /2 but have no well- 
defined charge. In D=2, fractionalization of the electron 
also occurs in the fractional quantum Hall effect |p7| , |l8| . 
For strongly-correlated systems, spin-charge separation 
has been proposed based on a slave-boson picture [|l9] p7| . 
However, recently Nayak |2j| has shown that the slave 
particles always remain confined. Our work suggests that 
rather than falling apart, electrons form clusters or com- 
posites when the Coulomb interaction is large and repul- 
sive. Such entities form pairs in the 2D Hubbard model. 
We propose that composite excitations of the kind we 
describe here offer a natural explanation for the absence 
of electron-like excitations in the ARPES experiments 
pTj on the underdoped cuprates. While it is tempting 
to draw a connection between pairing of composite ex- 
citations in the 2D Hubbard model and the pairing of 
composite fermions p9| in quantum Hall systems, this 
connection would be tenuous at best. 

Our approach then accomplishes two things: 1) First 
we are able to access the strong-coupling limit of the 
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Hubbard model and identify the relevant composite exci- 
tations. 2) We are able to show that such composites pair 
in the d x 2_ y 2 channel. Our work then represents a for- 
mulation of d-wave pairing from a microscopic model out 
of which emerges a hierarchy of new composite fermionic 
excitations. We find that the magnitude of the anoma- 
lous correlation function for composite excitation pairing 
increases as the compressibility in the normal state be- 
comes more negative. In the next section, we outline the 
essential details of the strong-coupling approach. In Sec. 
Ill, we derive the equations of motion for the compos- 
ite excitations and all the relevant correlation functions. 
The results for the anomalous pairing correlation func- 
tions are presented in Sec. IV. 



II. FORMALISM 

The starting point of our analysis is the on-site Hub- 
bard model 

H = - ]T Uj 4 a c aj + U n A n i ~ M Ul ( 4 ) 

i,j,<r i i 

where fy = t if are nearest- neighbour sites and 

zero otherwise, /i is the chemical potential and U the 
on-site Coulomb repulsion. Consider the retarded Green 
function 



((c ir ,c jl ))=6(t-t')({c ir (t),c jl (t')}) 



(5) 



where {A, B} is the anticommutator. As a result of the 
interaction term in Eq. (^), the equation of motion for 
this Green function will generate a new Green function 
that is proportional to {([c^ , H]; Cj i)) . This Green func- 
tion will contain the composite excitation r\i a . All higher- 
order Green functions will emerge from time derivatives 
of composite operators. This suggests that we should 
start with the Hubbard operators, rj a i = Ci a rii- a and 
£o-(i) = Ci a {\ — ni- a ) introduced previously. The util- 
ity of these operators is immediate. Their commutator 
with the interaction part of the Hubbard model yields 
— (ju — U)r\i a and — respectively. Hence, they can 
be used to diagonalise the t = Hubbard model. Conse- 
quently, if this basis is used to compute correlation func- 
tions, all higher-order correlation functions that will be 
generated will be multiplied by the hopping term. That 
is, these operators form the basis for a t/U <C 1 or equiv- 
alently a strong coupling expansion. The essence of this 
approach was first proposed by Linderberg and Ohrn ]3C| ] 
and has been applied to the Hubbard model J7 pl| , |3^ ] as 
well as the p-d model for the cuprates Q . 

To illustrate the utility of the composite operator tech- 
nique, we switch to the 4-component basis, 



/6t \ 

»7»T 



(6) 



Let us define 



m=iBtil>(i) = [fl>(i),ff\ 



(7) 



as the 'current' operator. Formally, we can write the n th 
element of the current as 



in(i) = Y K nm ip m (i) + 5j n (i). 



(8) 



We can project out the part of the correction, Sj, that is 
'orthogonal' to the composite operator basis by rewriting 



Sjn = J2{{Sjn^L})lJ^l+' 



(9) 



where 8<f> satisfies the equation, {Sip, = 0. Let us 
introduce the normalization matrix 



iu = ({>i>i,4}) = 



(2tt)2 

and the overlap matrix 



' 2 i- ik -( r i- r ! 



J(k) (10) 



d z ke 



M a (k) (11) 



where k-integration is over the Brillouin zone and Q the 
inverse area of the Brillouin zone. The elements of the 
I and K matrices contain only the mutual correlations 
among the constituents of the composite operator basis. 
Because {S(f>,^} — 0, we can write the overlap matrix 
as 



M = EI 



where 



E nm = K nm + ^({Sjn,^}) 1 ^- 



(12) 
(13) 



This matrix (the energy matrix) will play a crucial role 
in our approach. Because of the Sj contribution, the 
energy matrix contains higher order correlations which 
are not expressible as simple correlations of the ^-fields. 
It is from the Sj terms that the composite excitations 
arise. Consequently, the higher-order correlation func- 
tions must either be decoupled or used as self-consistent 
parameters in equations obtained by imposing various 
symmetries, such as the Pauli principle |^,^2|. As a re- 
sult of the orthogonality of 6<j> to the composite operator 
basis, we will neglect the contribution from 5cj> in all sub- 
sequent calculations. This approximation is equivalent 
to ignoring the dynamical corrections to the self-energy. 
Such contributions are expected to be small because the 
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Hubbard operator basis solves exactly the U = oo limit. 
The principal role of 5<p is to broaden the energy levels 
associated with the composite operator spectrum. The 
role of such dynamical corrections will be the focus of 
future study. Within this approximation, we can write 
any Green function 



in Fourier space as 

S(k,w) = (w-E(k))- 1 I(k) 



(14) 



(15) 



by using the equations of motion. This equation and 
the expression, M = EI, are the central equations of this 
approach. The M and I matrices involve correlation func- 
tions which can be obtained self-consistently and by in- 
troducing symmetry properties such as the Pauli princi- 
ple jj2j to close the equations of motion. It is the off- 
diagonal blocks of the energy matrix that contain infor- 
mation about anomalous pairing correlations. 

There are three types of self-consistent equations that 
enter this approach. The simplest of these involves the 
correlations between the ip fields that appear in M and 
I. The general expression for such correlation functions 
in terms of the corresponding Green function is given by 



<«)V>£(i)> = 



(2tt)3 



-1 



x I — J ImS mn (k, lo) 



(16) 



with f{u>) the Fermi-Dirac distribution function. The 
self-consistency of Eq. (jl^) follows from the dependence 
of E and I on correlations in the composite operator ba- 
sis. To obtain explicitly the self-consistent equations, we 
rewrite Eq. (|l5|) as 



5(k, W )= £ 

i=l,2 



where 



U> — 6j + irj U! + £j — it] 



A(±€i 



±2^ 



(17) 



(18) 



with chosen from (1, 2) but i ^ j, Cj the eigenvalues 
of the energy matrix and the matrix A is given by 



A(w) = DetH - E)H - E)~ X I. 



(19) 



The four eigenvalues correspond to the four composite 
operator bands: i=l,2 refer to the £ and t] bands, respec- 
tively, whereas +, — index the particle and hole states, 
respectively. If we now use Eq. ( |l7| ) in Eq. ( 16 ) , we obtain 
the general self-consistent equation, 



C m , n (i,l) = J d 2 ke^-^ [I mn + 

{( K t)->nn — (ft-i )mn) 



tanh ^ (20) 



for correlation functions within the composite operator 
basis where T = l/ksfi- 

It is well-known that in approximation schemes of this 
sort, certain correlation functions which vanish explicitly 
as a result of the Pauli principle, self-consistently iterate 
to a non-zero value. To avoid this shortcoming, we ex- 
plicitly maintained the symmetries [^],[52| imposed by the 
Pauli principle, namely 



(21) 



This is the second type of self-consistent integral equation 
that will be used in this method. Likewise, C2,i(i,i) — 
0. As a result of imposing these symmetry relations, 
we will find that our theory is completely particle-hole 
symmetric about half-filling. 

The final self-consistent equations arise from decou- 
pling correlations of composite fields that are not con- 
tained in the basis described by Eq. ([|). This procedure 
will be described in detail in the next section. 



III. COMPUTATIONAL MACHINERY 



A. Equations of Motion 



We start first by constructing the algebra 



1 



(22) 

(23) 
(24) 



of the composite operators, with for fj, = 1,2,3 are 
the Pauli matrices, a is the identity matrix, and 



5°""" = 



fJ-r, — 



nil 

~ C i1 C il 



(25) 



From these equations, it follows that the I-matrix is given 
by 



I = 



/1-f 




V A o 











f -Ao 
-AS 1-f 











A \ 
/ 



(26) 



where Ao — {c^Cii) and we have assumed that (rif) = 
{n±) = n/2. We now turn to the "current". We can 
construct the "current" from the equations of motion 
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j(i)=i^(i) = [iP(i),H] 
for the ^-fields with 

3 

J3(*)-/4i+E^ c ]l+M 



(27) 

(28) 
(29) 



and 



4f 



TTiJ. 



We also need the M-matrix. The distinct elements of this 
matrix are 



M 11 (k) = FT<{j 1 ,e}}> 
M 12 (k) 



-M(l - 2 ) - 4te 



-4to;(k)(l - n +p) 

,n 

2 



FT({j 1 ,r 7 ]})=4te-4ta(k)(--p) 

M 22 (k) = FT({i 2 ,r?]}) = -(a* - U)~ - 4te - Atap 

M 13 (k) = FT<{j'i,^}> = -4*7(k)fl + 8toA - 4tA c? (32) 

-( M + 4ta(k))A 
+4t 7 (k))0-4tA CI) 
(ji-U- 4to(k))A 
+4t 7 (k))6» + 4iA c? 
M 24 (k) = FT({j 2 ,»7j}> = -4i 7 (k)0 + 4iA c „ 



M 14 (k)= J FT({j 1 ,7 ? J) 
M 23 (k) = FT<{j 2 ,^}) 



with the normal correlations 

e = - (r/frjj), (33) 

and 

p = (n la n« a ) + (^(c^tH - ( CiTCii ( C t iC t T )«), (34) 
and anomalous correlations 



A 



and 



A cr) = (c iT r?f ; ) - (ci;r?f T ), 
= (citCij,[n T (ri +x) + ni(r, + x)]). 



(35) 



(36) 



(37) 



In these equations, FT signifies the Fourier transform as 
defined in Eq. (^), a represents an average over nearest- 
neighbour sites, x indexes nearest-neighbour sites, and 
a(k) = cos fez + cos k y . Unlike a(k) which is the coef- 
ficient of the normal correlation functions, the Fourier 



coefficient 7 (k) is a coefficient of an anomalous correla- 
tion function. Hence, it is sensitive to the sign change 
of the anomalous correlation functions as exp(ik • r) is 
summed over nearest neighbour sites. For s-wave sym- 
metry, a(k) — 7 (k) while for d-wave symmetry 7 (k) = 
cosfca, — cosfc y . The symmetry relationships among the 
elements of M are as follows: 



(38) 



The M13, Mi4, M 2 3, and M 24 contain the anomalous 
correlation functions in particular the 9 correlation dis- 
cussed in the introduction. 

With the M and I-matrices in hand, we now construct 
the energy matrix. From the structure of the normal- 
ization matrix and the symmetry properties of the M 
matrix (Eq. (|38|), it follows that the energy matrix is of 
the form, 



(30) 


Mil = 


-M 33 




Mis = 


-M 34 


(31) 


M 22 = 


-M 44 




M u = 


— M* 



E 



A 
B 



B 

-A 



(39) 



where A and B are 2x2 matrices. The off-diagonal 
blocks of the energy matrix determine the gap in the 
energy spectrum induced by pairing. To isolate the rel- 
evant parts of the normalization and M matrices that 
enter the energy gap, we consider the candidate pairing 
symmetries. No simplifications occur in the s-channel. 
However, we have verified that none of the anomalous 
correlation functions are non-zero either in this chan- 
nel. The same state of affairs occurs in p-wave symmetry. 
What about d-wave symmetry? Here several significant 
simplifications occur. First, as a result of the nodes in 
the gap, no purely on-site anomalous correlation func- 
tions survive. As a consequence, (c^cn) = Ao = 0. 
This leads to a vanishing of all the off-diagonal blocks 
of the normalization matrix. Consequently, only the off- 
diagonal blocks of the M-matrix enter the off-diagonal 
elements of the energy matrix. However, further simpli- 
fications occur. Consider for example, A c j and A cl) in 
Eqs. ([35j) and (|36|). These quantities are symmetrized 
and in addition involve a sum over the nearest neighbour 
sites in the Brillouin zone. Consequently, they vanish 
identically in the d-wave channel. The only anomalous 
correlation function that remains is 9, Eq. (|37|). At this 
level of theory, this is the only correlation function that 
does not vanish by symmetry conditions. This is the only 
anomalous correlation function that enters the gap in the 
energy spectrum. Hence, as advertised, any pairing will 
necessarily be governed by a non-traditional correlation 
function. This correlation function involves a composite 
operator that lives on a cluster of nearest-neighbour sites. 
By defining rji ]iy to be Ci a rij T , we rewrite Eq. (B7p as Eq. 
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(||) in which it is clear that pairing in 9 occurs between 
two composite excitations. 

B. Computational Procedure 

What remains to be done is the calculation of the 
anomalous correlation functions. Our assumption of sin- 
glet pairing implies that 

(c jT c a n 3T ) = (c iT c iJ .n J -_ r ) = (ct ; ct T n ir )*. (40) 

Consequently, we can define 9 as 

9 = 2(c n cnn jT ) (41) 

where the spin r is arbitrary. As this correlation function 
involves more than two basis operators, it must be decou- 
pled. We will follow a procedure analogous to that de- 
vised by Roth J31J in her treatment of the strong-coupling 
limit of the Hubbard model. To implement the Roth 
method, we consider the Green function 

H(i,j,t,t') = ((cJCi.fJjcr^OWy.O))- ( 42 ) 

Clearly, if this Green function is calculated, 9 can be 
obtained directly from an expression analogous to Eq. 
(|ilf). We proceed by constructing the series of Green 
functions, 

A a (i,j,k,t,t') = (<£T(M);c T 0\*0 t M*.O» 

B (y (i,j.M,*') = «»n(*.*);cT(j.*') t n«r(*,«')» ( 43 ) 
F a (i,j,k,t,i/) = ((£j(M);c T (j,t')V(M')}) 
G a (i,j,k,t,t') = ((? 7i t ^,<); CT 0-,<') t « CT (fc,t')}} 

in which the creation operator at time t is replaced by 
each of the composite operators in the 4-component ba- 
sis. The Green function of interest, H, is obtained by 
summing F and G and setting i = j. To calculate these 
quantities, we use the equations of motion for ip, Eq. (]?]), 
and in particular the approximation that j — idtip ~ Eip. 
Consequently, the equations of motion for the Green 
functions. 



id i 



( A a \ 




( A A 






F a 


= E(i) 


Ffj 


+ iS(t - £') 




\Gj 




\Gj 







(44) 



are directly related to the energy matrix in real space, 
E(i) and the term, f na = ({i/j n (i,t). l Cf(j,t)^n <T (k,t)}), 
that arises from the equal-time anticommutation of ip 
with the composite operators in the Green function. Note 
that in general, /„ is a linear combination of correlation 
functions in the composite operator basis and correla- 
tions associated with A + B or F + G. Consequently, we 



now have a closed set of equations from which 9 can be 
obtained. 

In the next step, we Fourier transform these equations 
of motion so that t — > u>, r& — r, — * ki, and Tj — Ti — > hi- 
Noting that uj — E = IS -1 , we obtain 



( A a (k 1 ,k 2l uj) ^ 
B (T (ki,k 2l uj) 
F a (k 1 ,k 2 ,uj) 

\ G a (ki,k 2 ,ui) J 



s(h + k 2l u)r 1 



/ hAh,k 2 )\ 

ha{k\,k 2 ) 



To extract 9 from these equations, we sum F and G, 
integrate over k 2 so that i = j and then using Eq. (|l^), 
we find that 



2( 



(45) 



where 



-4D 



n(2 — n) 

C = (Cn(i,j) + C 12 (i,j)) (C 13 (i,j) + C 14 (»,j)) 



(46) 



+- (C 22 (i,j) + C 12 (i,j)) (C 24 (i,j) + C u (i,j)) , (47) 



and D is the on-site double occupancy 
D = (n iT n a > = - - C 22 (i,i). 



(48) 



In the expression for 0, (i,j) denote nearest neighbour 
sites along the x-axis. We point out that Beenan and 
Edwards used a simple factorization scheme for 8. 
As these authors indicate 0, there is no unique way of 
performing this procedure. In our scheme, there are two 
distinct ways of decoupling the Green functions that lead 
to 9. That is, we could have started with another Green 
function 



H'(i,j,t,t,t') = << CT (M);ci(z,04aOc{(j,0>>- (49) 

Following the procedure outlined above, we arrive at the 
alternative expression for 9: 



9 = 



2C 



(50) 



with 

C' = : 



(C n (i,j) + C 12 (i,j)) (C 2i (i,j) + C u (i,j)) 



+- {C 22 (i,j) + C 12 (t,j)) (C 13 (i,j) + C u (i,j)) ■ (51) 
n 

In the following section, we compare the values of 9 ob- 
tained by both methods. 

Similarly it is possible to decouple the parameter p 
defined in Eq. ([mJ). Because Eq. ( p4[ ) contains only 
symmetric combinations of the fcrmionic operators, it 



G 



can be decoupled uniquely. From the procedure outlined 
above, we find that p is given by 



n 2 pt + <j>f% p x + p 2 p z 



1 + . 



(52) 



where 



Pi = 7^— (Cn(i,j) + C 12 (i, j)f + - (C 22 + C 12 f , (53) 
2 — n n 

P2 = 7^— ((C 13 (i,j) + C u (i, j)f + - (C 24 + Cuf (54) 



2- n 



(55) 



and 



A'? 



n(2 - n) 



(C n (i,j) + C 12 (i,j))(C 22 + C 12 ). (56) 



In the normal state, p 2 = as it contains only anomalous 
correlations and p reduces to the result derived by Roth 
pT| . The expressions for p and the self-consistent equa- 
tion for the correlation functions constitute the working 
equations of this method. 

To summarize, the equations for M and I contain four 
parameters, p, e, p and 9 that must be determined self- 
consistently. From Eq. ([53]) and the definition of the 
Hubbard operators, the self-consistent equations for e 
and n are 

e = C n (i,j)-C 22 (i,j) (57) 
n = 2(1 - Cu(i,i) - 2C 12 {i,i) - C 22 (i,i)). (58) 

In our computational procedure, we either used the two 
previous equations together with Eqs. ( [45| ) or (^0|). For 
the fourth self-consistent equation, we either used p or 
the equation that imposes the Pauli principle, Eq. (pl|). 



IV. RESULTS 

To set the stage for our results on the anomalous cor- 
relations, we discuss first the chemical potential. Shown 
in Fig. (Q) are three different calculations of the chemical 
potential as a function of the filling: 1) solid line- present 
method, 2) dashed-dotted line-Beenan and Edwards 0, 
and 3) dashed line-Avella and co-workers |3^| . Here, 
n = 1 corresponds to half-filling. The method of Avella 
and co-workers |52| is identical to the method used here. 
However, as is evident, our results are significantly dif- 
ferent. Particularly striking is the continuous increase of 
the chemical potential obtained by Avella and colleagues 
p2| as the filling is increased. This is in stark contrast 
our solution which levels off in the vicinity of half-filling. 
The difference between our treatments lies in that there 
are two solutions to the self-consistent equations. Avella 
and co-workers chose the solution that is higher in energy 



in the filling range 0.7 — 1.0. While the lower energy so- 
lution is the physically-relevant solution, this solution in 
the absence of pairing has a distinct negative compress- 
ibility in the vicinity of half-filling. It is for this reason 
that Avella and co-workers |32| criticized, the method 
used by Roth as it also gives rise to a negative compress- 
ibility as seen from the dashed-dotted line of Beenan and 
Edwards Avella and co-workers |32| advocated that 
imposing the Pauli principle in the self-consistent solu- 
tion to the equations of motion eliminates the negative 
compressibility. 




0.7 0.8 
Filling (n) 

FIG. 1. The chemical potential (measured in units of the 
hopping t) as a function of filling for U — 8t and T = 0. n=l 
corresponds to half filling. The dashed line corresponds to the 
work of Avella and colleagues |32| while the dashed -dotted 
line to the work of Beenan and Edwards j?| . Our work is the 
solid line and includes the effect of pairing. 



Our work shows that even if the Pauli principle is main- 
tained by means of Eq. (pl|), a second solution to the 
integral equations still exists along which the compress- 
ibility remains negative. To explore further the negative 
compressibility, we show in Fig. (^|) the role of pairing on 
the compressibility. The dashed line corresponds to the 
chemical potential in the absence of pairing while along 
the solid line 6* / 0. As is evident, pairing alleviates the 
negative compressibility almost entirely giving rise to a 
flattening of the chemical potential in the vicinity of half- 
filling as seen from the solid curve in Fig. (||). In fact, 
a key trend common to the curves shown in Fig. (Q) 
is that the compressibility is most positive when pairing 
vanishes. This result is particularly important because a 
negative compressibility occurs in numerous dilute elec- 
tron systems, such as the 2D electron gas for r s > 3 [[m). 
Our results also corroborate the earlier observation of 
Tandon, Wang, and Kotliar |}5| that the compressibility 
becomes negative in the U — > oo Hubbard model. For 
short-range Coulomb interactions, a negative compress- 
ibility signifies that the ground state of an electronic sys- 
tem is unstable relative to a uniform charge distribution. 
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Hence, a negative compressibility is typically associated 
with phase separation or stripe formation |36| ] . Our work 
explicitly shows that pairing alleviates this instability at 
least in the case of short-range Coulomb interactions. 
We have proposed that even in the case of long-range 
Coulomb interactions, pairing still obtains and alleviates 
the negative compressibility as well J37|]. 



0.04 



filling for the non-interacting system. This appears to be 
an accident of their approximations as our more accu- 
rate method shows that the peak in the order parameter 
occurs at a doping level of 10%. 




0.8 0.85 
Filling (n) 

FIG. 2. The chemical potential as a function of filling for 
U — 8t and T = in the absence of pairing (dashed line) and 
in the presence of pairing (solid line). The chemical potential 
is in units of t. The solid line clearly illustrates that pairing 
alleviates the negative compressibility. 



We now discuss explicitly the anomalous correlation 
functions. Shown in Fig. (J3J) are four different calcula- 
tions of the correlation function involving pairing of com- 
posite excitations. The two solid lines in Fig. (|J) were 
obtained from Eq. (^0|). On the upper curve, the Pauli 
principle was imposed whereas along the lower curve, the 
decoupling scheme was used to calculate the parameter 
p defined in Eq. (|34|). Hence, the lower curve corre- 
sponds to the method of Bccnan and Edwards 0. On 
the dashed curves, Eq. ( f45j ) was used to compute 9. Once 
again, along the upper dashed curve, the Pauli principle 
was used whereas along the lower curve, the decoupling 
scheme (for p) was used. Our results indicate that the 
anomalous correlations are largest and most stable when 
the Pauli principle is imposed by using the integral equa- 
tion, Eq. (pl|). When the decoupling scheme is used to 
obtain p, the two distinct decouplings for 9 yield vastly 
different results. This will result in huge fluctuations in 
T c as the work of Beenan and Edwards illustrates. 
What our work establishes is that if the Pauli principle 
is imposed in the self-consistent procedure outlined here, 
consistent pairing solutions exist regardless of the decou- 
pling scheme used to calculate 9. Note also that the two 
lower curves are peaked around a doping level of 20%. 
Beenan and Edwards associated great significance to 
this doping level as it corresponds to the filling at which 
the Fermi surface resembles the Fermi surface at half- 




1 1.2 1.4 

Filling (n) 

FIG. 3. Anomalous pairing correlation as a function of fill- 
ing for U = 8t and T = 0. The two solid lines were obtained 
from Eq. (^o|). On the upper curve, the Pauli principle was 
imposed whereas along the lower curve the decoupling scheme 
was used to calculate Eq. (|34|). On the dashed curves, Eq. 
( ^B| ) was used to compute 9. Once again, along the upper 
dashed curve, the Pauli principle was used whereas along the 
lower curve the decoupling scheme for p was used. Our results 
indicate that the anomalous correlations are largest and most 
stable when the Pauli principle is imposed. 



The remaining figures constitute the primary results 
of this method. Shown in Fig. (3|) is a calculation of 



9 in which the average of Eq. 45) and (|50|) as a func- 
tion of filling was used. Our results show clearly that 
d x 2_ y 2 pairing exists and such pairing is diminished as 
U increases. As our method is valid only in the large 
U limit, we cannot establish the lower value of U/t for 
which 9 is non-zero. Also of significance is the maximum 
in 9 at roughly 10% doping for the range of on-site re- 
pulsions studied. In addition, the peak in 9 moves to 
higher dopings as U increases. We emphasize that the 
pair-formation found here does not include the effect of 
phase fluctuations. Hence, the acutal doping regime over 
which true superconductivity with long-range order oc- 
curs might be significantly smaller than that obtained 
here. For example, if long-range order obtains at all, 
the optimum doping will appear shifted to higher doping 
values relative to the maxima of the curves shown here. 
Preliminary results by Manske, Dahm and Bennemann 
p8| using phenomenological spin-fluctuation approxima- 
tions offer some support of this conclusion. Also of note 
is the fact that our treatment preserves the particle-hole 
symmetry of 9. The temperature dependence of the order 
parameter is shown in Fig. (^). As is evident, the shape 
of 9 is characteristic of any order parameter that vanishes 
at a transition temperature. For U = 8i, we find that the 
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onset temperature, T* = .021i. In the cuprates t = .5eV. 
Hence, we obtain an onset temperature of T* = 100K for 
pair formation. Phase coherence occurs at a lower tem- 
perature, T c . As a consequence, the T* calculated here 
should serve as a realistic estimate of the psuedogap tem- 
perature |3£j . At optimum doping, T* should correspond 
to T c . Experimentally, optimal doping corresponds to 
roughly 20%. At this doping level and for U — 8t, we ob- 
tain that T c = 65K leading to a ratio of T*/T c = 1.5. For 
YBa 2 Cu 3 6 .95, T c = 92K while T* = HOif @. Typi- 
cally in the cuprates, 1.2 < T* /T c < 2 in the underdoped 
regime p9| . While our calculated values for T* and T c 
are not rigorous estimamtes, it is encouraging that they 
are not far off from the experimental values. 



0.01 




Filling (n) 

FIG. 4. Anomalous pairing correlation function as a func- 
tion of filling for various values of U for T = 0. We obtained 
taking the average of Eqs. (ph and @. 



V. FINAL REMARKS 

We have presented here an analysis ideally suited for 
the strong-coupling limit of the Hubbard model. Two 
key results are established. First, pairing in the Hubbard 
model occurs in the d x 2_ y 2 channel. Second, composite 
excitations rather than electon-like excitations pair to- 
gether. They live on a cluster of nearest-neighbour sites 
and are formed out of a hole and a singly-occupied site. 
Ultimately, such excitations could explain the absence of 
a well-defined electron-like peak in the ARPES experi- 
ments Jjrj] if the self-energy correction in Eq. ( |i~7| ) arising 
from the dynamical processes significantly broadens the 
energy levels near the Fermi energy. Sasaki, Matsumoto 
and Tachiki have evaluated such dynamical correc- 
tions in the context of the p-d model for the cuprates and 
found that the composite operator spectrum remains in- 
tact (thereby justifying the initial choice of the composite 
operator basis); however, broadening of all levels includ- 
ing those at the Fermi surface was observed. Matsumoto 
and Mancini |t2f have also performed similar calculations 



for the Hubbard model |42j and observed a broadening 
of the levels at the Fermi surface. 

0.01 1 > ■ ■ ■ 1 




0.005 0.01 0.015 0.02 0.025 0.03 
Temperature (T/t) 

FIG. 5. Temperature dependence of the anomalous corre- 
lation function at n = .9, the optimal doping as found in Fig. 
(4). All quantities are measured in units of t. The overall 
shape of 6(T) is consistent with the interpretation that 6 is a 
pairing gap. 

Longer-range Coulomb interactions, necessitate the 
retension of higher-order composite excitations. The sim- 
plest of such excitations will involve three sites. Comput- 
ing the equations of motion for the three-site composite 
excitations leads 4-site excitations. At each iteration of 
this procedure, the range of the composite excitations 
grows. We have verified that nearest-neighbour repulsive 
Coulomb interactions drastically diminish the anoma- 
lous pairing in 9. Trivially, attractive nearest-neighbour 
Coulomb interactions enhance pairing. However, to com- 
pletely settle the issue of pairing when repulsive next- 
nearest neighbour interactions are included, we must also 
retain the anomalous correlation functions that are gen- 
erated in the presence of the longer-range interaction. 
Our work here suggests that such correlation functions 
should be computed to determine if composite particle 
pairing survives in the extended Hubbard model. 

Nonetheless, within the on-site repulsive model, the 
correlation function (0) calculated here should be suffi- 
cient to describe the pairing process. Because 9 involves 
a product of the form f]ij\f]iji 1 the pairing mechanism is 
entirely local. From the form of 77^, it is tempting to 
conclude that pairing requires a doubly-occupied site to 
neighbour a singly-occupied site. In such a configura- 
tion, double occupancy can be shared between the two 
sites. However, this is just one of the many types of lo- 
cal configurations that gives rise to a non-zero value of 
9. If long-range phase coherence obtains, one can think 
of the condensate as a coherent superposition of all such 
resonating structures. This suggests that the pair-pair 
correlation function for the cexons should be calculated 
to verify one way or another if phase coherence obtains. 
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